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We present a method to solve the impulsive minimum fuel maneuver 
problem for a distributed set of spacecraft. We develop the method 
assuming a non-linear dynamics model and parameterize the prob- 
lem to allow the method to be applicable to multiple flight regimes 
including low-Earth orbits, highly-elliptic orbits (HEO), Lagrange 
point orbits, and interplanetary trajectories. Furthermore, the ap- 
proach is not limited by the inter-spacecraft separation distances and 
is applicable to both small formations as well as large constellations. 
Semianalytical derivatives are derived for the changes in the total AV 
with respect to changes in the independent variables. We also apply 
a set of constraints to ensure that the fuel expenditure is equalized 
over the spacecraft in formation. We conclude with several exam- 
ples and present optimal maneuver sequences for both a HEO and 
libration point formation. 

Introduction 

In this paper, we present a direct approach to solve the impulsive, minimum fuel ma- 
neuver problem for a distributed set of spacecraft. To equalize the fuel expenditure among 
spacecraft, we enforce a set of nonlinear constraints. We present an explanation of the 
method and the mathematical theory assuming a general, nonlinear dynamics model that 
can be expressed in an inertial or rotating coordinate system depending on the problem 
being solved. This ensures the method is applicable to a wide range of spaceflight regimes. 
Furthermore, since a nonlinear dynamics model is used, the method is not limited to small 
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inter-spacecraft separations and is equally applicable to large constellations and close for- 
mations. 

Trajectory optimization has a rich history and many techniques have been developed 
over the last fifty years. Most techniques can be classified as either direct or indirect. 
Betts 1 presents an excellent survey of trajectory optimization techniques. Guzman 2 et.al 
present a survey of indirect methods. The approach developed here is a direct method 
and is an extension of techniques that have their origins in research performed during 
the Apollo era 3,4 and later extended by D’Amario 5 et. al . , and Howell 6,7 et. al. Their 
are several new contributions contained in this work. First, the method is generalized to 
permit minimum fuel optimization for a maneuver sequence involving a set of m spacecraft. 
We have also generalized the method to find the optimal launch vehicle injection orbit to 
minimize fuel during the initial spacecraft deployment phase. Finally, we have reformulated 
the cost function to remove a naturally occurring singularity in the gradient, without loss 
of generality. 

We begin this paper with a mathematical problem statement defining the cost and 
constraint functions. The cost function is the sum of the total AV of all spacecraft over 
an entire maneuver sequence. The parameterization of the problem is discussed for two 
types of maneuver sequences: initialization and reconfiguration. Next, we discuss how to 
evaluate the cost and constraints. The approach requires solving a number of Initial Value 
Problems (IVP) and Two Point Boundary Value Problems (TPBVP). However, given the 
speed of modern computers, the method is surprisingly fast. Next, the gradient of the cost 
function and the Jacobian of the constraints are derived. We discuss numerical issues in 
implementing this approach. We conclude the paper with several test problems in different 
flight regimes to demonstrate the applicability of the method. A brief explanation of the 
notation is included in Appendix 3. 

Problem Statement and Parameterization 

There are numerous ways to pose the minimum fuel formation maneuver problem. We 
assume that the desired relative motion is driven by mission requirements and has been 
determined a-priori. The goal is to develop a technique to achieve the desired relative mo- 
tion in a minimum fuel manner. We define the minimum fuel problem for a formation of m 
spacecraft, where the trajectory of the k th spacecraft has n k total maneuvers, as 

Given: r OA; , v 0fc , t Qk , r/ fc , v /fc , t fk 

where r Gfc , v 0fc , and t Qk are the initial position, velocity, and epoch of the k th spacecraft 
respectively, r f k and v f k are the final position and velocity respectively of the k th spacecraft 
at a reference epoch t f k 

Solve: 

min(J) (1) 

where 

m 

( 2 ) 

k = 1 
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where m is the number of spacecraft and AV& is the total weighted AV expended by the 
k th spacecraft and is calculated using 

Tlk 

AV k = Y,f°( Av ^ Av i k (3) 

3 = 1 

where A Vjk is the j th maneuver on the k th trajectory. The scalar function f s (Avjk) in Eq. (3) 
is included to remove a singularity in the derivative of J for small A Vjk- This function will 
be discussed in detail in a later section. For now it suffices to say that f s (Avjk) = 1 for 
values of A Vjk large enough not to cause numerical problems. When appropriate, we impose 
the following set of constraints to equalize the AV among the spacecraft. 


(AVi - AV 2) 2 

< 

toll 

(4) 

{AV 2 - av 3 ) 2 

< 

tol2 

(5) 

(AV^-i - AV m ) 2 

< 

tolm—\ 

(6) 

(AVi - AV m ) 2 

< 

tolui 

(7) 


where toll is the desired tolerance for constraint one and so on. 

We assume a dynamics model given by the second order differential equation 

r = f (r, r, t) (8) 

The formulation of the cost and constraints, and their derivatives, is performed for the 
general nonlinear dynamics model seen in Eq. (8). We discuss the specific dynamics model 
used for software implementation and validation in a later section. 

At the heart of solving an optimization problem is the problem parameterization. In 
general, the cost function J shown in Eq. (2) is a scalar function of vectors that we can 
write as 

m 

J = Y, AVk = J ( X ’ C) (9) 

k = 1 

where X is the vector of independent variables being manipulated by an optimization rou- 
tine, and C is a vector of constants associated with the problem. The constraints in Eqs. (4) 
- (7) can be written as 

G = diag (ZAVAV t Z t ) = G(X, C) < TOL (10) 

where AV = [A Vi AV 2 ... AV m _i AV m ] T , Z is constant matrix discussed in a later sec- 
tion, and TOL is a vector of tolerances associated with the constraints. Before attempting 
to solve the problem we must choose which variables and degrees of freedom associated 
with the problem are to be included in X, which will be varied by a numerical optimization 
routine, and which degrees of freedom are to be included in C and treated as constants. 
The choice of X and C will influence the convergence properties of the numerical routine 
and determine the types of problems the method can solve. 
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In the remainder of this section, we discuss two parameterizations to solve the impulsive 
minimum fuel maneuver problem for distributed spacecraft. For convenience, we categorize 
the types of maneuver sequences as initialization sequences or maintenance/reconfiguration 
sequences. In an initialization sequence, all spacecraft begin at a common time and location 
on the same orbit. An initialization sequence is performed to take each member spacecraft 
from the common parking orbit to its desired final location defined by a final orbit and time. 
In a maintenance maneuver sequence, the spacecraft are not initially located on the same 
orbit and hence are not collocated before the maneuvers are performed. The fuel optimal 
maintenance sequence is a more general problem and we begin by discussing the approach 
taken here. 


Maintenance and Reconfiguration Sequence 


An illustration of a formation reconfiguration sequence is shown in Figure 1. For sim- 
plicity the diagram only illustrates a two-spacecraft sequence. Before the beginning of the 
maneuver sequence, each spacecraft is located on a unique orbit. We define the initial orbit 
for the k th spacecraft, V ok , as the locus of points that are the solutions to the initial value 
problem 

**1 k = (j'lk^ok-) [^ofc) Vofc]) (11) 

for varying values of t\ k , the time at which the k th spacecraft performs its first maneuver. 
Recall that r ok and v ok are the position and velocity respectively, of the k th spacecraft at 
some reference epoch t ok . Throughout the paper we use 0 r to denote the position portion of 
the solution to the Initial Value Problem (I VP) of the differential equation shown in Eq. (8). 
The initial conditions to the I VP shown in Eq. (11) are contained in square brackets. The 
initial time is the second argument and the final time is the first argument. 


The goal of the reconfiguration sequence is to depart from the trajectories defined by 
V Q k , and move each spacecraft to a new location on Vfk which is defined as the locus of 
points that are the solution to the initial value problem 

** nick = 0 0'ni c kitfki {^fk-t^fk}) ( 1 ^) 

for varying values of t nkk , where rik is the number of maneuvers on the k th trajectory, and 
Yfk and Vfk are the desired final position and velocity respectively, of the k th spacecraft 
at some reference final epoch tf k • According to this definition, t nkk is the time of the last 
maneuver of the k th spacecraft. 


Between the initial and final maneuvers, the k th spacecraft can perform intermediate 
maneuvers at locations rj k given by 


I 2 < j <n k - 1 
\ 1 < k < m 


(13) 


Note that the reason 2 < j < n k - 1 is because and r nkk are determined from the initial 
value problems shown in Eq. (11) and (12). Because we allow the time of the first and last 
maneuver for each spacecraft to vary, as well as the times of the interior maneuvers, we see 
that ti kl the maneuver times, are given by 


f 1 < i <n k 
\ 1 < k < m 


(14) 
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There may be numerical difficulties if is not ordered such that t 2 k > tik • This is discussed 
in a later section. 

From Eqs. (11-14) we see that the variables that define the maneuver sequence are: tik > 
t 0 k , v 0 k,tfk-> Y fk, and tik . We can now parameterize the problem by choosing 

which variables to include in X and which to include in C. For the general minimum fuel 
reconfiguration problem we choose to include as independent variables the time of the first 
maneuver of the k th spacecraft, t\k, the times of the interior maneuvers, t *•*., and the times 
of the final maneuvers, t nk k • We also include the locations of the interior maneuvers, r jk, 
as independent variables. We treat the state components t Q k , r 0 k , and v 0 *. that define the 
initial orbit, and the state components tfk , r/fc, and Vfk that define the final orbits, as 
constants. In summary, for the maintenance and reconfiguration problem we can write 

^ f 2 < j < n k - 1 


x = [tik r Jkf 

< 1 < i 

< n k 

(15) 

{ l < k 

< m 


C=[t ok v T ok 

T T 

Vok tfk r /fc 

T 1 T 

V A J 

(16) 


The parameterization of the reconfiguration problem shown in Eqs. (15) and (16) is 
chosen such that the boundary conditions at rx^ and r nk k are satisfied implicitly. The 
optimizer is not tasked with satisfying another set of complicated non-linear constraints to 
satisfy the boundary conditions along the initial and final orbits V Q k and Vfk respectively. 
The down-side to this approach is that for each cost function evaluation we must solve a 
number of initial value problems (IVP) and two point boundary value problems (TPBVP). 
Solution of the TPBVPs is discussed in a later section. In the next section we discuss the 
parameterization of the initialization sequence. 
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Initialization Sequence 

The second type of maneuver sequence we address is an initialization sequence. In an 
initialization sequence, all spacecraft initially lie on a common orbit defined by the locus of 
points V 0 , as shown in Figure 2. Hence, for the initialization sequence we have 

*ok = r o ? v ok = v 0 , t 0 ^ — t Q and , V Q k V Q (17) 

We refer to V 0 as the injection orbit and it may or may not be known a-priori. For large 
formations that may require significant fuel expenditure from all spacecraft in formation 
to achieve the desired relative geometry, it may be desirable to include the injection orbit 
itself as an independent variable. In this case, the solution to the optimization problem 
yields the minimum fuel initialization trajectory for each spacecraft as well as the optimal 
launch vehicle injection orbit. Referring to Figure 2, we see that although all spacecraft 
initially lie on V Q , it is not required that all spacecraft perform their first maneuver in the 
initialization sequence at the same epoch and hence the spacecraft can depart from V Q at 
different locations. The location of the first maneuver for the k th spacecraft is determined 
by the epoch of its first maneuver t\k, and the state vector that defines V 0 • We assume that 
V 0 is unique trajectory defined by the locus of points that is the solution to the following 
I VP for varying values of t\k\ 

rik = ko) ^o]) ( 18 ) 

where t Q , r G , and v 0 , are the launch vehicle injection conditions, and t\k is the time of the 
first maneuver of the k th spacecraft. There is a unique final trajectory for the k th spacecraft 
defined by the locus of points, Vfk, which is the solution to the following IVP for varying 
values of t nk k : 

r nfc/c = [ r //e? v //c]) (^®) 

We assume the values tfk, rjk, and vy* have been determined a-priori to maximize the 
mission return. Hence these values are constants for the initialization problem. 


6 


For complete generality, we assume that the epoch of the initial maneuver for each 
spacecraft and the epoch of the final maneuver to insert the spacecraft into its location 
in formation can freely vary. Note that the initial epoch of the first maneuver of the k th 
spacecraft, t\k is not the same as t Q . The epoch t Q is part of the state vector that uniquely 
determines the injection orbit, while t\k is the epoch when the k th spacecraft departs from 
the injection orbit. We also assume that the times and locations of the interior maneuvers 
for all spacecraft can freely vary. For the initialization problem we define X and C as 
follows: 


Uk rJ k ] T | 

[ii 

3 

i 

i—l 

1 

£ g 
VI VI 

(20) 

l i< 

k 

< m 


C=[t fk rj k \ 

r T ] T 
r fk \ 



(21) 


The parameterization of the cost and constraint functions in an optimization problem 
greatly influence the convergence properties of the method. In this section, we discussed 
a parameterization of the minimum fuel maneuver problem for two types of maneuver se- 
quences for a set of distributed spacecraft. Given the problem parameterizations developed 
in this section, we are ready to look into the details of evaluating the cost and constraint 
functions. 

Cost and Constraint Function Evaluation 

In this section we take the problem parameterization given by X and C in the last 
section and discuss how to calculate the cost J(X, C) and the constraints G(X, C). Calcu- 
lating the cost and constraints is similar for both the initialization and the reconfiguration 
problem and so we address them simultaneously. The first step is to redimensionalize the 
vector of independent variables X. Once we have a dimensionalized set of variables given by 
X and C, we solve a set of I VPs to obtain the position vectors that define the locations of 
the first and last maneuvers of each spacecraft. The next step is to solve a set of TPBVPs 
that yield A Vj^. Finally, we solve for the cost and constraints using Eqs. (2), (3) and (32). 
Each step is discussed in more detail below. 


Step 1: Redimensionalize X 

It is often necessary to work with nondimensional variables when solving optimization 
problems using numerical methods. We begin by assuming that the numerical optimiza- 
tion process uses non-dimensional variables where the transformation from dimensional to 
nondimensional variables is given by 

X' = diag(X/) -1 X — X s (22) 

where X' is the nondimensional form of X, diag(Xj) is a matrix with X/ on the diagonal 
and zeros for the off-diagonal terms, and X j and X A . are vectors of the same length as X. 
We assume that diag(Xy) -1 exists. The inverse transformation is simply 

X = diag(X / )(X' + X 5 ) (23) 


7 


It is more convenient to work in dimensional variables to evaluate the cost and constraints 
and the redimensionalization is performed before step 2 is performed. X f and X 6 . are chosen 
according to the numerical issues associated with the particular problem being solved. 

Step 2: Solve IVPs 

To evaluate the cost and constraints, we need to be able to calculate 

A^H|v+*-v7J { \ ^ l ^ (24) 

where a superscript denotes a condition immediately before an impulse maneuver, and 
a superscript “+” denotes a condition immediately after an impulsive maneuver. The 
quantities and are determined by solving the following initial value problems 

v u = <t> v (hk,tok, [r 0 fc, v ofc ]) (25) 


V nfcfc — ^ tfki [ r /fci v /fc]) 

where <fi v denotes the velocity portion of the solution to the I VP, and tfk , r/fc, and are 

contained in C and treated as constants. If we are solving a maintenance problem, r 0 fc, 

and v 0 /~ are also contained in C. However, if we are solving an initialization problem then 
the initial boundary conditions may be contained in X. To solve the TPBVPs found in the 
next step we will need iq k and r nfc fc. However, these are simply the position portions of the 
solution to the IVPs in Eq. (25) and (26). 

Step 3: Solve TPBVPs 

The remaining components of the impulsive maneuver vectors needed to calculate Eqs. (24) 
are determined by solving a set of TPBVPs. From Step 2 we found and r nAs jt, and the 
variables r ^ and tjk are contained in X. Therefore, we know the positions and times for all 
maneuvers and we have N Lambert’s problems to solve where N is given by 

771 

N = ^(n fe - 1) (27) 

k= 1 

There are numerous well-known approaches to solving Lambert’s problem. We use a simple 
multiple shooting method outlined below. The algorithm is described as follows where we 
drop the subscript fc, for now, to simplify the notation. Given an initial position and an 
initial velocity v^-, both at time tj , find 8vj applied at time tj so that we achieve rj+i at 
tj- |_i. Figure 3 illustrates the problem. The dark black arc denotes the path a spacecraft 
would follow if no 8vj were applied. For this arc, the final position denoted by r a is not 
equal to the desired final position Tj+i. Hence, for the dark black are, 5rj+i ^ 0. The 
dashed arc denotes the trajectory that is the solution to Lambert’s problem. For this arc 
rj+i = r a , or £rj+i = 0. To solve for 8vj such that 8rj+i = 0 first define Xj as 
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Figure 3 Illustration of One Trajectory Arc 


then <5x^+1 is given by 


Sx j+ i = $ (t j+ 1 , tj) 5xj 


where tj + 1 and tj are fixed. We can write (tj+\,tj) as 


^ {tj+litj) 


Aj+i,j Bj+ij 


(29) 


(30) 


We can solve for Av^ such that is zero by iterating on Eq. (31) until (r J+ i — r a ) meets 
a user defined tolerance. 

Avj = B jl hj (r j+ i - r a ) (31) 

Upon convergence, for each trajectory arc we save vj~ , v“, and $ for use in calculating the 
cost, constraints, and their derivatives. It is important to note that this approach assumes 
that B j exists. We address cases when this is not true in a later section. Using the above 
algorithm, we solve all N Lambert’s problems. Knowing the solution to each trajectory arc 
permits the calculation of all maneuvers A Vjk and with this information we can evaluate J 
and G. 


Step 4: Solve for J and G 

From the previous steps we know v^. and v~ k . We use Eqs. (24) to calculate A Vjk- Next 
we can evaluate f s (Avjk) using Eqs. (41)-(43). J can be calculated using Eqs. (2) and (3). 

The constraints in Eqs. (7) are more conveniently written for mathematical manipulation 
as 

G = diag (ZAVAV t Z t ) < TOL (32) 

where TOL is a vector of tolerances, AV is given by 

AV = [A Vi AK 2 ... AV m ] T (33) 
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and Z is given by 


/-I 1 0 0 ... 0\ 

0-11 0 ... 0 

0 0 -1 1 ... 0 


0 ... -1 1 0 0 
0 ... 0 -1 1 0 

0 ... 0 0-11 

\1 0 0 ... 0 -1 / 


The above equations allow us to calculate J(X, C), and G(X, C). To take full advantage 
of the power of a numerical optimization routine, it is also useful to calculate the gradient 
of J and the Jacobian of G. We devote the next section to this topic. 


Cost and Constraint Derivatives 


Gradient-based numerical optimization routines such as Sequential Quadratic Program- 
ming (SQP) perform best when supplied with analytic derivatives of the cost and con- 
straints. Providing derivatives is often nontrivial, yet it is advantageous because it increases 
the speed of convergence by requiring fewer function calls and avoids numerical problems 
associated with finite differencing. 


Below we derive the gradient of the cost function and the Jacobian of the constraint 
functions. We need the derivatives of the cost and constraint functions with respect to 
the independent variables shown in Eqs. (15) and (20). Specifically, we require analytic 
expressions for 


dJ_ 

dx 


d_ 

dx 


i >n 


k=l 


(34) 


and 


<9G „ 

— = 2 diag 


^ZAV 


<9AV t 

dx 



where x is a generic independent variable and 


(35) 


dAV _ [AVi dAV 2 
dx dx dx 


dAVm 


dx 


(36) 


In Eqs. (34)- (36) we see terms of the form dAV k /dx which can be written 


5^ -££/.(**>** 

j= 1 

One can show that Eq. (37) can be written as 

dAVk _ / df s (Avjk) f s {Avjk)\ ^Aw- k 

dx dAv jk + Av jk J dx Jk 


(37) 


(38) 
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By inspecting Eq. (38), we see that we must determine dAvJ k /dx. This term is non- 
trivial and in the next few subsections we derive this term for the specific independent 
variables used in this work. The term f s {Avj k )/Avjk can cause difficulty for small Avj k . 
Traditionally, to determine the total AV, we would set f s = 1 and we have the new 
relationship 

AV{ = Y,^v jk (39) 

3 = 1 

However, if f s = 1 then we have the following term that appears in the derivative of AV k . 

fs(AVjk) = 1 ( 4 Q) 

A Vjk Avjk 

This term is obviously singular for small Avj k . The difficulty this causes is significant 
because the optimization process tends to make Avjk as small as possible and sometimes 
Avjk can approach zero. We can mitigate problems caused by small values of Avjk by 
carefully choosing f s (Avj k ). 

By inspecting the term f s (Avj k )/Avj k , we see that we need f s (Avj k ) to approach zero 
as fast or faster than Avj k so that the limit as Avj k — > 0 is defined. There are many possible 
choices for such a function. We have chosen to partition the function f s (Avj k ) into three 
regions. For small values of A Vj k , where Avj k < A vl we define f s {Avj k ) as 

fs{ Av jk ) = Av 2 jk (A Vj k < A v L ) (41) 

For intermediate values of Avjk, where Avl < Avjk < Avh, we define f s (Avjk) as 

f s (Avjk) = aAvj k + bAvj k + cAv 2 jk + dAv jk ( Av L < Av < Av H ) (42) 

For large values of Avj k , where Avh < A Vj k , we define f s (Avj k ) as 

f s (Avj k ) = 1 {Avh < Av jk ) (43) 

The values of Avl and Avh m the above equations are chosen depending upon the numerics 
of the particular problem being solved. However, Avl must always be less than Avh- 


The form of f s described above removes the singularity in the derivative of f s {Avj k ) for 
small values of Avj k . As Avj k — > 0 we see that where Avj k < Avl 

to to = o (44) 

^Vjk Avj k >0 Avj k 

and we have a non-singular function for all values of Avj k . The four constants in the quartic 
function in Eq. (42) are chosen so that f s (Avj k ) and its first derivative are continuous. The 
derivation of the constants a, 6, c, and d is shown in Appendix 1 and the expressions are 
found in Eqs. (114-117). 


The second term in Eq. (38) that is non- trivial is 


dAv Jk 

dx 


( 45 ) 
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where z is a dummy variable that represents an arbitrary component of X given in Eqs. (15) 
and (20). First let’s identify some derivatives that are zero after looking at the physics of 
the problem. From inspection of Figure (1) and (2) we see that 


dAvf k 

— — — = 0 if n 7 ^ k 
or (n 

(46) 

®AvJ k 

3 =0 if 

(47) 


where n is dummy index used to denote the trajectory number. Eqs. (46) and (47) come 
from the fact that changing the time or position of a maneuver on one trajectory does not 
affect the AV along another trajectory. For this reason we can drop the subscripts n and 
k without loss of generality. Hence the remaining derivatives to be calculated are 

9AvJ dAvJ 
dvg ’ dti 

where the second subscript on each variable is assumed to be k. For convenience, we break 
down the derivatives shown above into three categories and devote a subsection to each. 
The first category contains initial boundary derivatives which are derivatives with respect 
to £ 0 , r 0 , v 0 , and t\. The second category contains derivatives that are with respect to 
internal maneuver times and locations, tg (tg € ti) and rg respectively where 2 < i < — 1 . 

The third category contains derivatives with respect to the final time, t Uk . The next three 
subsections discuss the three types of derivatives in detail. 


Initial Boundary Derivatives 

We define the initial boundary derivatives to be those with respect to t 0 , r 0 , and v D . 
By inspection of Figure ( 1 ) we can see that 


SAv£ 
dt i 


dAvJ dAvJ 


dt 0 


dr 0 


dAvJ 

= 0 if j > 2 

<9v 0 


(48) 


This is due to the fact that changing ro does not change AV 3 . Likewise changing time t Q 
does not change AV 3 . Therefore the nonzero initial boundary derivatives are 


9Av^ dAv2 dAvJ 3Av^ <9Av^ SAv^ <9Avj" SAv^ 

dt\ 5 dt\ ’ dt 0 7 dt 0 5 dr Q ’ dr 0 ’ dv Q ’ dv 0 


Let’s begin by looking at the derivatives with respect to t\. Assume we have the TPBVP 
illustrated in Figure 4 and defined by 


Given: t 0 , r 0 , v G , t 2 , r 2 
Find: vf such that r \ t2 = r 2 
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Figure 4 Illustration of Initial Boundary 


We know the solution has the following form 

n = 4 > r [r2, V2 ]) 


(49) 


v+ = <p v [r 2 ,v 2 ]) 


(50) 


The symbol <j> denotes the solution to the initial value problem with initial conditions 
contained in square brackets. The second time argument is the initial time and the first 
time argument is the final time. So, in Eq. (49), T 2 and v^" are the initial conditions with 
an initial time of and a final time of t\. Hence, in this case we are back-propagating 
because £2 > £ 1 . The superscript r and v denote the position and velocity portions of the 
initial value problem respectively. By inspection of Figure 4, we see that changing t\ will 
result in a change in ri, vj~, v+, and v^\ Taking the derivative of ri and with respect 
to t\ yields 


dr\ d(p r d4> T dv 2 

dt\ dt\ dV dt\ 

(51) 

dv+ _ d<p v d<p v dv 2 
dti dti * dV dt\ 

(52) 


where V is a dummy variable to denote differentiation with respect to velocity. These 
equations can be rewritten as 


V 1 =v^ +B tl ,t 2 


<9v 2 
dt 1 

dv f _ a + , D ^ 2 " 
— - ai +D tl>t2 — 


(53) 


(54) 


Solving the system of equations, and using Bqs. (118) and (119) from Appendix 2, we obtain 

dv~ 


2 _ 


dt i 


= (Ct 2ltl ~ U t2>tl B t2 * tl At 2; tj) (v, - v+) 


dvf , 

1- = a+- R 

dt i 


(v, -v+) 


(55) 

(56) 
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By inspection, we see that 

<9v7 

dh _&1 

(57) 


(58) 

The derivatives with respect to r 0 and v G can be derived by starting with 

a solution of 

the form 


r 2 = <j> T (t 2 ,ti,[ri,vf}) 

(59) 

v 2 = 4> v {h,h, [ri, vjj"]) 

(60) 

where we note that changes in r Q or v G result in changes in iq, vj", vf, and 
derivatives of Eqs. (59) and (60) with respect to r 0 yields 

. Taking 

dr2 d(j) r dr\ d(/) r dv+ 

dr Q dR dr Q 1 dV dr Q 

(61) 

d^2 d<P v & r i ^ v i" 

dr Q dR 0r o 1 dV dr Q 

(62) 

where R is dummy variable to denote differentiation with respect to position, 
tions can be rewritten as 

dwt 

0 = A t2,fl A tl,t° + Q r 

These equa- 

(63) 

aV 2 C \ i D dV ' 

(64) 

Finally, solving for the desired derivatives yields 


Q r ^*2 >*1 "^2 >*1-^1 i*o 

(65) 


(66) 

Taking derivatives of Eqs. (59) and (60) with respect to v 0 yields 


ch *2 d<p r dr\ dcj) r <9v+ 

dw Q dR dv Q * dV dw Q 

(67) 

dw^ d<p v dri d(j) v dv+ 
dv Q ~ dR dv 0 + dV dv Q 

(68) 

These equations can be rewritten as 


dvt 

0 - A t2,tl B «l,to + B t 2 ,tl 

(69) 
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( 70 ) 


dr 0 

solving for the desired quantities yields 


dv 7 dvt 


dr Q 


d y i = _ R -i A R 


(71) 


dv. 


dv„ 


2 _ 


= ( C t2,tl “ D t2,tl B i 2 !tl A *2,tl) 


(72) 


Finally, the last two nonzero derivatives with respect to r c and v„ can be written simply as 

<9v7 


0r o 

dv7 


L = c tl|i 


to 


d\ 0 


= D 


t\,to 


(73) 

(74) 


The last initial boundary derivatives to be determined are the derivatives with respect 
to t Q . The derivatives of ri and vj“ with respect to t Q can be found by starting with a 
solution of the form 

r 0 = <f> r (t 0 ,ti, [ri,vf]) (75) 


v 0 = 4> v [ri,v 1 ]) (76) 

where we note that changes in t a result in changes in ri, vj", v+, and • Taking the 
derivative of Eqs. (75) and (76) with respect to t Q yields 


dr Q d(f) r d<fi r dr i ^ <90 r dw x 


dt 0 dt Q dR dt Q dV dt Q 


dv Q d(f> v d(j) v dr\ d(j) v <9v : 

dt Q dt Q ^ dR dt Q dV dt Q 


These equations simplify to 


. . dri dv 1 

o =v » +A< --9(: +B « i atr 

_ dr i dv7 

0 = a o + C Wl — + D to , fl — 


solving for the desired derivatives yields 

dr i 


of = ~ ( A h,to V o + B <i,to a o) 


dv 

dt, 


— (^-'ti,t o y 0 d" B ti,t 0 a o) 


(77) 

(78) 

(79) 

(80) 

(81) 

(82) 
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Table 1 Initial Boundary Derivatives with Respect to t Q 

Term 

d t 0 


dvr 

dvZ 

dv+ 

~ Dt 2 ,t 1 D tl) t 2 ) _1 Ct 2 ^ 1 (A ili t 0 v 0 B £l)io a G ) 
— (I - D^ 2j f 1 D( lj t 2 ) _1 Ct 2iil (At 1 ^ 0 v 0 + B^^a^) 

0 


Recall that changing t Q changes ri, v 1 , v+, and v 2 . We can find the 
dv+ /dt Q and dw 2 /dt 0 starting with solutions of the form 

expressions for 


v+ = <p v {tut 2 ,[r 2 ,v 2 ]) 

(83) 

V 2 = 0 t, (^2,ti,[ri,vj t ’]) 
Taking the derivative with respect to t Q yields 

(84) 


dv+ d<fi v dr 2 d<p v dw 2 

dt Q dR dt Q * dV dt Q 

(85) 


dv 2 d(\> v 3 yi d(p> v <9v+ 

dt Q dR dt 0 * dV dt 0 

(86) 


This can be rewritten as the system of equations below 


dv+ 

dt a 


- D tl , t2 


dv 2 

dt a 


(87) 


dvo 

Hu 




dvf 

dt 0 


Solving the system of equations for the desired quantities yields 


dt 0 


(I _ D( 2itl D fl)t2 ) 1 C t2) { 1 (A tlil<) v 0 + B tli t 0 a 0 ) 


( 88 ) 


(89) 


= -(I - D f2)fl D tl , t2 )- 1 C t2 , tl (A tl , to v 0 + B tl ,< 0 a 0 ) (90) 

Note that Eqs. (89) and (90) contain an element of the inverse of the state transition matrix. 
See Appendix 2 for a discussion of how this can be calculated without inverting the entire 
6 by 6 state transition matrix. 

We have completed the derivation of the initial boundary derivatives and they are sum- 
marized in Tables 1-3. Now we move on to the internal derivatives. 
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Table 2 Initial Boundary Derivatives with Respect to t\ 


Term 

dt 1 

3vf 

a r 

<9v+ 

a l" - (vf - v+) 

<9V2 

( C t 2 ,h ~ D «2,tl B ( 2 |ti A t2,tl) ( V 1 - V l") 

dv+ 

0 


Table 3 Initial Boundary Derivatives with respect to r Q and v Q 


Term 

dr Q 

dv a 

avf 

Cti,to 


dyf 


— Ai 2.tl B <l,to 

dvz 

( C <2,il “ B t^t, A -t2,<l ) A tl,to 

( C t2,tl - D f2,tl B ^|ti A <2,«l) B tl,to 

dv+ 

0 

0 


Internal Derivatives 


The internal derivatives are defined as derivatives of Avj with respect to the internal 
maneuver times and positions, or: 


and 


8 A vj 1 

r i< 

j 

< n k 

(91) 

dr e I 

l 2< 

t 

<n k -l 

dAv J 1 

r i< 

j 

< n k 

(92) 

dte 1 

l 2< 

l 

<n k - 1 


By inspection of Figure (1) we can see that some of the derivatives shown in expressions 
(91) and (92) are zero. For example, changing r3 does not change Avi or Avs. Likewise 
we see that changing the time ^3 does not change Avi or AV5. In general we can write 


<9AvJ 

_ 3 =0 if j < i - 1 or j > £ + 1 
or? 


(93) 


<9AvJ 

. ■ 3 =0 if j < £ - 1 or j > t + 1 
dte 


(94) 


In other words, changing the time or position of the j th maneuver only changes the following 
maneuvers: Avj_i, Avj, and Av^+i. Thus, the only non-zero internal derivatives are 


dAvJ_ x dAvJ dAvJ +1 
dre ’ dre ’ dve 

The non-zero internal derivatives can all be determined by closely investigating Lambert’s 
problem. Recall that Lambert’s problem is to find vj and v^, given 17, tj, rp and tp. 
We see that the non-zero internal derivatives can be determined from the derivatives of the 
solution to Lambert’s problem , vj and v^, with respect to changes in 17, tj, rp and tp. 
These derivatives appear in the literature, 5,7 and they can also be found by the approach 
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Table 4 Time Derivatives of the Solution to Lambert’s Problem 
dti dtp 


dvf a+ + B J it/ Aj Fit ,v+ _ 

dv]r (DtF,t/B tF 1 ,t r A tF ,t J - Q F , tj )v/ a F - Dt^Bj^Vj, 


Table 5 Position Derivatives of the Solution to Lambert’s Problem 


dr j 

dr p 

dv^ 

dvp 

Ct F ,ti ~~ 

■p> — l 
a t Fi ti 


used in this paper to determine the derivatives of cost and constraints with respect to the 
independent variables. Table 4 contains a summary of the derivatives of and with 
respect to changes in the initial and final time tj and tp respectively. Table 5 contains a 
summary of the derivatives of and with respect to changes in the initial and final 
position r/ and r p respectively. The left hand column of the tables contain the numerator of 
the derivative definition and the horizontal titles contain the denominator of the derivative 
definition for a particular derivative. So, for example, we see that 


<9v+ 

dti 


= »/ + 


(95) 


where and A tF)tj come from the STM of the trajectory that is the solution to Lam- 

bert’s problem, v* is the velocity immediately after the first impulse, and is the accel- 
eration immediately after the first impulse. Using the derivatives in Tables 4 and 5, we can 
determine the remaining internal derivatives which are summarized in Table 6. 


It is interesting to note that the derivatives dvj /dte and /dte are explicit functions 
of the accelerations and a^T. For spacecraft flying in environments with nonconservative 
forces these terms may be significant. However, when we evaluate Eq. (38), the accelerations 
appear in pairs such as (a^ - - aj). This term is identically zero when the flight regime 
consists only of conservative forces and so acceleration terms do not come into play unless 
there are forces that are explicit functions of the spacecraft velocity. 


Final Boundary Derivatives 

The final boundary derivatives are defined as derivatives of AvJ with respect to t nk . 
From inspection, we see that 


dAvJ 

dt nk 


= 0 


if j < n — 1 


(96) 
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Table 6 Interior Derivatives 


Term 


dte 


dr e 

dv e-i 


0 


0 

dy ti 


-B 


■D — 1 

dv e 



-l V « 


dv+ 

a e 


,te v e 


dv e+i 



C te+ute)vt 


dv+i 


0 


0 


Table 7 Final Boundary Derivatives with Respect to t nk 


Term 

dtn k 

d Vn - 1 

0 

dvt-i 

1 ( V «fc _ V «fc) 


a n k + (v+ - V" ) 


< 


Following a similar approach as in the last two sections, the final boundary derivatives are 



(97) 

d <-i _ R -i / + v - x 

D n fc ,nfc- 1 Wn k n k ) 

UL n k 

(98) 

9w n _ + 

dt nk nk 

(99) 

~ a nfc ~^~^ri ky n k -l^n k ,n k -l ( V n/t V n k ) 

UL n k 

(100) 

The quantity v+ fc is found using 


v n fc =0” (*»*.*/>[ r /> V /]) 

(101) 


and and come from Eq. (8) evaluated at the pre and post maneuver conditions 
for maneuver n^. The derivatives with respect to t nk are summarized in Table 7. This 
completes the derivation of the required derivatives. We now move on to discuss numerical 
issues and implementation. 

Numerical Issues and Implementation 

Subtle issues in the implementation of a numerical optimization approach often dra- 
matically influence the speed of convergence and the quality of the solutions. The method 
presented here has several numerical issues that can be accommodated if they are under- 
stood by the analyst. In this section, we address these numerical issues as well as discuss 
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details of the software implementation of the method used to solve the test cases in the next 
section. The first difficulty occurs when certain components of the STM are not invertible. 
The second difficulty is an artifact of the problem parameterization and occurs when the 
times of two maneuvers are the same. The third difficulty is when there are small AVs 
along a trajectory. 

The most obvious numerical difficulty appears in calculating the derivatives of the cost 
and constraints. The equations for the derivatives in the previous section assume that 
the components of the STM are invertible. Stern 10 showed three cases in which certain 
components of the STM are not invertible for the two-body problem. They are: 

1. The difference between the initial and final times is an integer multiple of the orbit 
period. 

2. The difference between the initial and final true anomalies is given by Ntt, for N = 
0,1,2, 3, . . . (Note that case 1 is a subset of case 2) 

3. A certain trigonometric function of the eccentricity and the eccentric anomalies is 
zero. 


The first two cases cause surprisingly less difficulty than one might initially expect. Test 
cases on the Hohmann transfer, whose solution involves a transfer angle of exactly 7r, have 
shown that components of the STM are invertible for transfer angles within le-4 degrees of 
the singularity that occurs at the solution. The result is that the method finds solutions 
within le-6 m/s of the known analytic solution. Furthermore, this method has been devel- 
oped to solve real-world problems where perturbations are included in the dynamics model. 
In these cases, the solution rarely occurs exactly at transfer angles (2N + l)n. Finally, the 
third case occurs for multiple revolution solutions. The examples investigated in this work 
did not contain multiple revolution trajectory arcs. 

Another subtle singularity arises from the parameterization of the problem. The opti- 
mization routine controls the independent variables shown in Eq. (15) for a maintenance 
or reconfiguration sequence and the independent variables shown in Eq. (20) for an ini- 
tialization sequence. Suppose we are optimizing a maneuver sequence that contains two 
spacecraft (m = 2) and that each spacecraft performs four maneuvers ( n\ = n 2 = 4). For 
this scenario, £22 and £ 32 , the times of the second and third maneuver of the second space- 
craft respectively, and r 2 2 and r 32 , the location of the second and third maneuvers of the 
second spacecraft respectively, are varied freely by the optimizer. It is possible that during 
the optimization process, that the optimization routine may try the following: 

£22 = £32 r 22 7^ r 32 (102) 

There is no solution to Lambert’s problem for the above conditions. Furthermore, as £ 22 — > 
t 32 , AV 22 — ► 00 . There are several ways to handle this difficulty. One is to apply a set of 
constraints such that > c where c is a constant chosen according to the problem 

being solved. However, if during the optimization process this problem occurs, it is likely 
that one of the maneuvers is not necessary. The simplest solution is to stop the optimization 
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process, remove one of the maneuvers, in this case either the maneuver at time £22 or 1 32, 
and restart the optimization process. 

The third numerical difficulty is the singularity that occurs when f s {Avjk) = 1 and 
A Vjk = 0 and was discussed in previous sections. By defining f s (Avjk ) as in Eq. (41) - (43), 
we can remove this singularity. Physically, the form chosen for f s (Avjk) removes small 
values of Avjk from the cost function. Care must be taken in defining the constants Avl 
and Avh according to the problem being solved. 

The formulation of the cost function, constraints, and gradients was performed without 
regard to the specific dynamics model chosen. However, in a software implementation, 
we must choose a dynamics model including the reference frame in which to express the 
equations of motion, and the forces and perturbations to include. We have chosen to work 
in the Earth Mean J2000 Equatorial system. We have chosen to include accelerations from 
the spherical Earth, J2, and third body point mass accelerations from the Sun and Moon. 
The resulting dynamics are given by 


^21 = ~G 


(m 2 + mi) 


r 21 


' 21 


,.3 \ ’n 


m. 


~ r j2 + a J 2 




(103) 


where a subscript “1” represents the spacecraft, a subscript “2” represents Earth, a subscript 
“3” represents the sun, a subscript “4” represents the moon, and r2i is the position of the 
spacecraft with respect to the Earth and so on. The term a j 2 is given by 


a J 2 


/ 5xz 2 \ 


3»/2/ie-^e 
2 r 5 


X — 


y- 


3z 


r 2 

r*2 


(104) 


where x, y, and 2 are the inertial components of the spacecraft position vector, and fi e and 
R e are the Earth’s gravitational parameter and radius respectively. The STM is calculated 
by defining 

x = [ r T v T ] 7 (105) 

then 



(106) 


then in first order or state form 



(107) 


The derivatives in Eq. (107), while non-trivial, are well known and not presented here. 
Finally, the differential equation governing the STM is given by 


4> = A$ 


(108) 


The differential equations shown in Eqs. (103) and (108) are numerically integrated as a 
coupled system of 42 first-order ordinary differential equations. 
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The remaining issues involved in the software implementation are concerned with opti- 
mizer selection, convergence criteria, scaling, and bookkeeping of the considerable amount 
of data required to evaluate the cost function and its gradient and the constraint func- 
tions and their Jacobian. We have chosen to work in MATLAB and therefore have used 
the fmincon SQP algorithm available in the MATLAB Optimization Toolbox. For speed, 
all numerical integration of the equations of motion and the STMs is performed by using 
MATLAB mex functions compiled from C code. Likewise, all TPBVPs are solved in mexed 
C functions. The TolX, TolFun, and TolCon convergence tolerances of fmincon were set to 
le-12 for most problems and we refer the reader to the MATLAB Optimization Toolbox 8 
documentation for more information on the definition of these settings. 

The final issue involves bookkeeping. For each cost function evaluation, we require an 
initial guess in order to start the TPBVP solver discussed in a previous section. For the first 
function evaluation, the code uses an initial guess provided by the user. For subsequent cost 
function evaluations, the code uses solutions from previous cost function evaluations as the 
initial guess. This assures that the initial guesses for the TPBVPs evolve with the optimizer 
iterations. There is an STM matrix associated with each segment of each trajectory and 
these STMs are calculated when the TPBVPs are solved. The STMs for each segment are 
saved, and then used in the calculation of the Gradient and Jacobian. Great care must be 
taken to save the STMs in a convenient manner and use them correctly in the equations for 
the Gradient and Jacobian. 

We now move on to discuss several example applications. 

Applications and Examples 

In recent years, numerous distributed spacecraft missions have been proposed in a diverse 
set of flight regimes and employing a wide range of inter-spacecraft separation distances. 
The examples we have chosen demonstrate the applicability of the method to different 
flight regimes and different types of distributed spacecraft missions. In example 1, case 1, 
we choose not to solve for the optimal injection orbit. In example 1, case 2, we do solve 
for the optimal injection orbit. Similarly, in some cases we have chosen to apply the AV 
equalization constraints, and in other cases we have chosen not to apply the constraints. 

The first example is a Highly Elliptic Orbit (HEO) formation of four spacecraft that 
forms an ideal tetrahedron with sides of 100 km at apogee. For this HEO formation, we 
apply the method to initialize the tetrahedron in the presence of perturbations for three 
cases: 


1. Minimize J, do not vary injection orbit, do not apply AV equalization constraints. 

2. Minimize J, vary the injection orbit, do not apply AV equalization constraints. 

3. Minimize J, vary the injection orbit, apply AV equalization constraints. 

The second example is a formation in a large-amplitude Lissajous orbit about the Sun- 
Earth L 2 point. The formation is composed of two spacecraft that initially have the same 
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state and epoch. The final configuration is a separation of 200 km between the two space- 
craft. For this example, we present an optimal initialization sequence to demonstrate the 
applicability of the method to multi-body flight regimes. 

Example 1: HEO Formation 

Recall that example 1 is a HEO formation of four spacecraft that forms an ideal tetra- 
hedron of 100 km at apogee. The desired final orbit states are found in Table 8. The initial 
orbit, called the injection orbit, is found in column five. The desired final states for the 
spacecraft are found in columns one through four. The injection orbit was determined by 
back-propagating the desired final states for approximately 24 hours, and averaging the 
states of the four spacecraft. A three maneuver initial guess sequence was then created by 
solving the resulting four lambert problems and adding a small intermediate AV for a total 
of three maneuvers on each of the four trajectories. The AV s for the initial guess found 
using this approach are shown in Table (9). We see that the total AV for the initial guess 
is 256.63 m/s. 


Table 8 States Defining Desired Orbit for Example 1 


Property 

S/C 1+ 

S/C 2t 

S/C 3 T 

S/C 4* 

Injection Orbit 1 

a (km) 

42095.74 

42095.74 

42095.74 

42095.74 

42095.70 

e 

.8181803 

.8181822 

.8162243 

0.8169300 

0.8173807 

i (deg.) 

18.00000 

18.00000 

18.02006 

17.94851 

17.99215 

u) (deg.) 

90.00012 

90.09161 

90.04603 

90.04597 

90.04581 

(deg.) 

0.1020429e-3 

0.1026448e-3 

.1024955e-3 

0.1023461e-3 

0 

v (deg.) 

178.6875 

178.6709 

178.6699 

178.6733 

179.9916 


t 22 Mar 2012 12:00:00.0000 f 23 Mar 2012 11:17:34.2939 


Table 9 AV’s for Initial Guess for Example 1, Case One 


Property 

S/C 1 

S/C 2 

S/C 3 

S/C 4 

Avi (m/s) 

24.78 

21.87 

33.98 

35.01 

Av 2 (m/s) 

6.97 

7.058 

5.964 

6.357 

Avz (m/s) 

20.04 

22.44 

36.87 

35.28 

Y2 Av i ( m A) 

51.79 

51.37 

76.82 

76.65 


Total AV = 256.63 m/s 


For example 1, case 1, an optimal maneuver sequence was found for the HEO forma- 
tion by selecting the independent variables X and the constants C as shown in equations 
Eqs. (15) and (16). Hence, we treated this as a reconfiguration problem where initially 
the spacecraft were collocated at the same point on the initial orbit. We did not vary the 
state of the initial orbit, and we did not apply the AV equalization constraints. The AV s 
associated with the optimal solution for case 1 are found in Table 10. We see the total AV 
is 71.476 m/s and maximum accumulated AV for a single spacecraft is about 20 m/s. This 
is a fuel savings of 185 m/s over the initial guess. 
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Table 10 Optimal Solution: Example 1, Case 1 


Property 

S/C 1 

S/C 2 

S/C 3 

S/C 4 

A^i (m/s) 

8.864 

8.783 

5.372 

6.593 

Av2 (ra/s) 

1.475 

2.258 

3.276 

6.33 

Avs (m/s) 

5.797 

6.438 

8.899 

7.572 

v j (m/s) 

15.956 

17.479 

17.5466 

20.495 


Total AV = 71.476 m/s 


For case 2, we began with the converged solution to case 1 as the initial guess, but 
allowed the initial orbit to be included in the independent variables as opposed to being 
considered a constant. Hence, X and C were defined as shown in equations Eqs. (20) 
and ( 21 ). The solution yields a minimum fuel maneuver sequence, and the optimal launch 
injection orbit. The AVs for case 2 are shown in Table 11 and the optimal launch injection 
orbit in Table 12. The total AV for case 2 is 20.79 m/s. We see that varying the injection 
orbit has a dramatic influence on the total AV. There is a 51 m/s improvement over case 
1 and the maximum accumulated AV for a single spacecraft is about 7.3 m/s. 


Table 11 Optimal Solution: Example 1, Case 2 


Property 

S/C 1 

S/C 2 

S/C 3 

S/C 4 

A^i (m/s) 

0.247 

4.187 

1.003 

4.5e-5 

Av 2 (rn/s) 

0.523 

0.593 

2.382 

7.212 

A 7/3 (m/s) 

3.836 

4.1e-5 

0.791 

3.6e-5 

I2 Av j (m/s) 

4.606 

4.780 

4.176 

7.212 


Total AV = 20.774 m/s 


Table 12 Example 1, Case 2 Optimal Injection Orbit 


Property 

Injection Orbit^ 

a (km) 

42091.69 

e 

0.8168908 

i (deg.) 

17.999768 

u (deg.) 

89.78643 

fi (deg.) 

0.181559 

v (deg.) 

180.14348 


t 22 Mar 2012 11:58:35.4160 


For case 3, we choose the same independent variables as example two, however we applied 
the AFequalization constraints shown in Eq. (32). We used the solution from Example 1, 
case 2 as the initial guess for case 3. The AVs for an optimal solution for case 3 are shown 
in Table 13. Notice that the total AFfor each spacecraft is 6.093 m/s. Comparing to case 
2 we see that enforcing the AKequalization constraints resulted in a small penalty of about 
3.34 m/s. We now move on to discuss Example 2 . 
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Table 13 Optimal Solution: Example 1, Case 3 


Property 

S/C 1 

S/C 2 

S/C 3 

S/C 4 

Avi (m/s) 

0.789 

5.151 

1.460 

0.002 

Av 2 (m/s) 

2.460 

1.940 

3.254 

6.068 

Av^ (m/s) 

2.844 

0.002 

1.379 

0.023 

E Au i ( m / s ) 

6.093 

6.093 

6.093 

6.029 


Total AV = 24.116 m/s 


Example 2: Libration Point Formation 

The second example consists of a two-spacecraft libration point formation problem about 
the Sun-Earth L 2 point. The desired states of the spacecraft after the maneuver sequence, 
and the injection orbit and initial epoch are shown in Table 14. The states are for a large- 
amplitude Lissajous orbit. The initial guess maneuver sequence for Example 2 contains 


Table 14 States Defining Desired Orbit for Example 2 


Property 

S/C 1 T 

S/C 21 

Injection Orbit* 

x (km) 

76863.56215 

76915.15497 

-70150.29946 

V (km) 

-407755.1147 

-407729.4652 

1358773.158 

z (km) 

-33798.44889 

-33989.96985 

983153.4317 

x (km/s) 

-0.04802618 

-0.048026184 

-0.10552622 

y (km/s) 

1.12799544 

1.127995438 

-0.00980635 

z (km/s) 

0.13812967 

0.138129667 

-0.05853318 


t 11- Jan-2004 09:35:33.6099 J 09-Jul-2004 09:35:33.6099 


four maneuvers for each of the two spacecraft. The two internal maneuvers were computed 
by first finding a two-maneuver trajectory and then placing two small maneuvers, spaced 
equally in time, along the original two-maneuver trajectory. The AVs associated with the 
initial guess for Example 2 are found in Table 15. The total AV for the initial guess for 
each spacecraft is about 2.6 m/s. 

We chose to treat this problem as a formation maintenance problem and we did not vary 
the initial orbit. Furthermore, due to the fact that the AV of the initial guess is already 
small, we did not apply the AV" equalization constraints. 

The AV^s for the optimal maneuver sequence for Example 2 are shown in Table 16. 
Comparing the optimal solution to the initial guess we see that the optimal solution is 
about a one order-of-magnitude improvement over the initial guess. Also, while we did not 
enforce the Ay equalization constraints, the Ays for the optimal solution are nearly equal 
at approximately 0.109 m/s. 
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Table 15 Ay's for Initial Guess for Test Ca se Two 


Property 

S/C 1 

S/C 2 

Av\ {m/s) 

0.2272 

0.2272 

Av 2 {m/s) 

1.171 

1.171 

Av$ {m/s) 

1.098 

1.087 

Av± {m/s) 

0.1766 

0.28 

E Av j 

) 2.6728 

2.7651 

Total AV 

= 5.4379 m/s 

Table 16 Optimal Solution: Example 2 

Property 

S/C 1 

S/C 2 

Av\ {m/s) 

0.05679 

0.0568 

Av 2 {m/s) 

0.04973 

0.04972 

Avs {m/s) 4.058e-005 

4.066e-005 

Avi {m/s) 

0.00279 

.00275 

E Av j ( m / s ) 

0.10935 

0.10932 

Total AV 

= 0.21867 m/s 


Conclusions 

Distributed spacecraft and formation flying missions have been proposed for numerous 
flight regimes and for a wide range of interspacecraft separations. It is useful for the mission 
analyst to have at his or her disposal a method that is applicable to as many mission 
scenarios and dynamics regimes as possible to perform optimal maneuver planning. 

In this work we presented a direct approach to find minimum fuel maneuver sequences 
for distributed spacecraft missions. The cost function is defined as the cumulative AV of 
all spacecraft in formation and we proposed a set of optional constraints to equalize the fuel 
expenditure among spacecraft over a particular maneuver sequence. The method requires 
solving a set of IVPs and TPBVPs for each cost function evaluation. However, given the 
speed of modern computers the method is not prohibitively slow. Analytic derivatives of 
the cost and constraints were derived to take full advantage of the power of the numerical 
optimization routine. The method was applied to two test problems: a HEO formation and 
a libration point formation. Several optimal scenarios were presented. 

Their are several new contributions to the literature contained in this work. We gener- 
alized methods previously developed in Refs [3] - [7], to permit minimum fuel optimization 
for a set of m spacecraft. We also generalized the method to find the optimal launch ve- 
hicle injection orbit to minimize fuel during the initial spacecraft deployment phase. The 
cost function was reformulated to remove a naturally occurring singularity in the gradi- 
ent, without loss of generality. We also formulated a set of constraints to equalize the 
fuel expenditure among spacecraft. These modifications, together with the work performed 
by previous researchers, provides an optimization technique for minimum fuel distributed 
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spacecraft maneuvers in multiple flight regimes including LEO, HEO, libration and in- 
terplanetary trajectories. Furthermore, the method is not limited to small interspaceraft 
separations and is applicable to small formations or large constellations. 
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Appendix 1 


In Table 17, we see the three forms that the function f s (Avj/) can assume depending 
on the magnitude of A vjk. The constants in the quaxtic equation are chosen so that f s and 
its first derivative are continuous for all values of A Vjk- If we label the different portions 
of f 3 as Fi, F 2 , and F3 as shown in Figure 5, then the conditions to ensure a continuous 
function and its first derivative are 


dFi 


dAv 


Av l 


dF 2 


dAv 


A VL 


(109) 
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Table 17 Definition of f s (Avjk) for different values of A Vjk 

Avjk 

fs(Avjk) 

A v jk < Avl 
Av l < Avj k < Av h 
Av h < Av^ 

aAVj k + bAvj k + cAv 2 k + dAvjk 
1 


Fi(Avl) — F 2 (Avl) 


(110) 


dFj, 

dAv 


— 0 

Avh 


(in) 


F-2 {Avh ) = 1 

These conditions yield the following system of linear, algebraic equations 


[4 Avl 

3A y ^ 

2Av l 

i N 


(a\ 


(2Av l \ 

Avl 

A»i 

Av 2 l 

Avl 


b 


Avl 

4Av 3 h 

3Av% 

2A vh 

1 


c 


0 

\A v 4 h 

Av 3 h 

Av 2 h 

Av h ) 

w 


V i / 


The solution to this system of equations is 

AvlAv 2 h + A vl + Av^ - 3 Avh 
° (-A v H + Avl) 3 Av 2 h 


( 112 ) 


(113) 


(114) 


Av 2 l Av 2 h + Av\ + Av^Avl - 2AvlAvh + Avjj — 2Av 2 H 
(-Avh + Avl) 3 Av 2 h 


(115) 


Av\ + Av\A vh + ^Av 3 H Av 2 L + AvjfAvL - &AvlAv 2 h + Av b H 
(-Av h + Avl) z Av 2 h 

d = 2 (A vl + Avjj - 2 Avh)Av 2 l 
Avh(~Avh + Avl ) 3 


(116) 

(117) 


Appendix 2 


In order to calculate all of the required partial derivatives, it is sometimes necessary 
to calculate the inverse of the state transition matrix. Using well known formulas for the 
inverse of a block matrix, and assuming that all of the necessary inverses exist, one can 
show that 


tj) 


Ai+i,j Fj+ij 

Cj+ ij F)j+ij _ 


(Aj+ij Bj + ijDj +1 jCj+ij) 
(Bj+i,j ~ -A-j+ijCj+ijDj+ij) 


(Cj+i,j F)j B- +l] A j + ij) 

(Bj+ ij — Cj +ijA.j + ijBj+ij) 
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From the above, we determine that 


D 


B 7d + i 

B 7l 


jj+l-B j t j+l 


C j + 1 , j Dj+ijBj+ij Aj + 1 j 
~ B j+i,j A j+ 1 d 


(118) 

(119) 


The relations above can be useful in reconciling the various expressions for the partial 
derivatives found in the literature. It should also be noted that for Hamiltonian systems, 
the state transition matrix is symplectic 9 which implies that 


*(*i+i»*i) 1 = 




c 


j+hj 


j+ 1,J 
■®i + lj 


-l 


n T 

—C* T 

lj 




Appendix 3: Notation 

r Position vector 

v Velocity vector 

m Number of spacecraft in formation 

rik Number of maneuvers along k th trajectory 

State transition matrix 
A Upper left 3x3 partition of 

B Upper right 3x3 partition of 

C Lower left 3x3 partition of 

D Lower right 3x3 partition of <3> 

Av^ j th impulsive maneuver on k th trajectory 
A vjk Magnitude of j th maneuver on k th trajectory 

V 0 k Initial trajectory of k th spacecraft 

V/k Final trajectory of k th spacecraft 

X Vector of independent variables 

C Vector of constants 

i Internal maneuver location index, 2 < i < rik - 1 

j Maneuver time index, 1 < j < rik 

k Trajectory index, 1 < k < m 
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